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O . ABSTRACT 
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^ ! Using an analytical model and numerical simulations, we show that acoustic 

Q ! waves generated by turbulent motion in intracluster medium effectively heat the 

i central region of a so-called "cooling flow" cluster. We assume that the turbulence 

^ i is generated by substructure motion in a cluster or cluster mergers. Our analytical 

^ i model can reproduce observed density and temperature profiles of a few clusters. 

O i We also show that waves can transfer more energy from the outer region of a 

I cluster than thermal conduction alone. Numerical simulations generally support 

O ! the results of the analytical study. 
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1. Introduction 



^ ' For many years, it was thought that radiative losses via X-ray emission in clusters of 

(73 I galaxies leads to a substantial gas inflow, which was called a "cooling flow" (Fabian 1994, 
and references therein). However, X-ray spectra taken with ASCA and XMM-Newton fail 
to show line emission from ions having intermediate or low temperatures, implying that the 
cooling rate is at least 5 or 10 times less than that previously assumed (e.g. Ikebe et al. 
1997; Makishima et al. 2001; Peterson et al. 2001; Tamura et al. 2001; Kaastra et al. 2001; 
Matsushita et al. 2002). Chandra observations have also confirmed the small cooling rates 
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(e.g. McNamara et al. 2000; Johnstone et al. 2002; Ettori et al. 2002; Blanton, Sarazin, & 
McNamara 2003). 

These observations suggest that a gas inflow is prevented by some heat sources that 
balance the radiative losses. There are two popular ideas about the heating sources. One 
is energy injection from a central AGN of a cluster (Tucker & Rosner 1983; Bohringer & 

Morfin 1988; Rephacli 1987; Binney & Tabor 1995; Sokcr ct al. 2001; Ciotti & Ostriker 2001; 
Bohringer et al. 2002; Churazov et al. 2002; Soker, Blanton, & Sarazin 2002; Reynolds, Heinz, 
& Begelman 2002; Kaiser & Binney 2003). Recent Chandra observations show that AGNs at 
cluster centers actually disturb the intracluster medium (ICM) around them (Fabian et al. 
2000; McNamara et al. 2000; Blanton et al. 2001; McNamara et al. 2001; Mazzotta et al. 2002; 
Pujita et al. 2002; Johnstone et al. 2002; Kempner, Sarazin, & Ricker 2002), although some 
of them were already discovered by ROSAT (Bohringer et al. 1993; Huang & Sarazin 1998). 
Numerical simulations suggest that buoyant bubbles created by the AGNs mix and heat the 
ambient ICM to some extent (Churazov et al. 2001; Quilis, Bower, & Balogh 2001; Saxton, 
Sutherland, k Bicknell 2001; Briiggen k Kaiser 2002; Basson k Alexander 2003). The other 
possible heat source is thermal conduction from the hot outer layers of clusters (Takahara k 
Takahara 1979, 1981; Tucker k Rosner 1983; Friaca 1986; Gaetz 1989; Bohringer k Fabian 
1989; Sparks 1992; Saito k Shigeyama 1999; Narayan k Medvedev 2001). 

However, it has already been known that the ICM heating by AGNs or thermal con- 
duction has problems. For the AGN heating, the efficiency of the heating must be quite 
high (Fabian, Voigt, k Morris 2002). Moreover, the intermittent activity of an AGN makes 
the temperature profile of the host cluster irregular, which is inconsistent with observations 
(Brighenti k Mathews 2003). For the thermal conduction, stability is the most serious prob- 
lem; either the observed temperature gradient disappears or the conduction has a negligible 
effect relative to radiative cooling (Bregman k David 1988; Brighenti k Mathews 2003; Soker 
2003). Moreover, thermal conduction alone cannot sufficiently heat the central regions of 
some clusters (Voigt et al. 2002; Zakamska k Narayan 2003). Although a "double heat- 
ing model" that incorporates the efi^ects of simultaneous heating by both the central AGN 
and thermal conduction may alleviate the stability problem (Ruszkowski k Begelman 2002), 
Brighenti k Mathews (2003) indicates that the conductivity must still be about 0.35 ± 0.10 
times the Spitzer value. 

In this paper, we consider another natural heating source. In the ICM, fluid turbulence 
is generated by substructure motion or cluster mergers. From numerical simulations, Nagai, 
Kravtsov, k Kosowsky (2003) showed that the turbulent velocities in the ICM is about 20%- 
30% of the sound speed even when a cluster is relatively relaxed. Such turbulence generates 
acoustic waves in the ICM. Compressive characters of the acoustic waves with a relatively 
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large amplitude inevitably lead to the the steepening of the wave fronts to form shocks. As 
a result, the waves can heat the surrounding ICM through the shock dissipation. A similar 
heating mechanism has also been proposed in the solar corona; the waves are excited by 
granule motions of surface convection (Osterbrock 1961; Ulmschneider 1971; McWhirter, 
Thoncmann, & Wilson 1975). The idea of wave heating in the ICM was proposed by Pringle 
(1989), but the study was limited to order-of- magnitude estimates. In this paper, we study 
the wave heating by an analytical model and numerical simulations. We use cosmological 
parameters of Qq — 0.3, A = 0.7, and Hq — 70 km s~^ Mpc~^ unless otherwise mentioned. 



2. Analytical Approach 
2.1. Models 

In the ICM the magnetic pressure is generally negligible against the gas pressure (Sarazin 
1986). Therefore, acoustic waves (strictly speaking, fast mode waves in high-/? plasma) could 
carry much larger amount of energy than other modes of magnetohydrodynamical waves. 
We expect that turbulence in the ICM excites acoustic waves that propagate in various 
directions. In this paper, we focus on the acoustic waves traveling inward, which play an 
important role in the heating of the cluster center. These waves, having a relatively large 
but finite amplitude, eventually form shocks to shape sawtooth waves (N-waves) and directly 
heat the surrounding ICM by dissipation of their wave energy. We adopt the heating model 
for the solar corona based on the weak shock theory (Suzuki 2002; Stein & Schwartz 1972). In 
this section, we assume that a cluster is spherically symmetric and stationary. The equation 
of continuity is 

M = -Anr^fw , (1) 

where M is the mass accretion rate, r is the distance from the cluster center, p is the gas 
density, and v is the gas velocity. The equation of momentum conservation is 
dv ^ GM{t) _ldp 1 1 d 

dr r2 pd.r pc^l + [(7 + l)/2]a^} r2 ^ ^' 

where G is the gravitational constant, M(r) is the mass within radius r, p is the gas pres- 
sure, Cs is the sound velocity, 7(= 5/3) is the adiabatic constant, and is the wave velocity 
amplitude normalized by the ambient sound velocity (a^ = Sv^/cg). For the actual cal- 
culations, we ignore the term vdv/dr because the velocity is much smaller than the sound 
velocity except for the very central region of a cluster where the weak shock approximation 
is not valid {ayj > 1; see §2.2). The wave energy flux, F^, is given by 

^ —pclal U + ^^-^ayA . (3) 
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Note that the sign of equation (3) is the opposite of equation (7) of Suzuki (2002), because we 
consider waves propagating inwards contrary to those in Suzuki (2002). The energy equation 
is written as 
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where ks is the Boltzmann constant, T is the gas temperature, //(= 0.61) is the mean 
molecular weight, tjih is the hydrogen mass, is the electron number density, and A is the 
cooling function. The term V ■ indicates the heating by the dissipation of the waves. We 
adopt the classical form of the conductive flux for ionized gas. 



dr 



(5) 



with kq = 5 X 10^^ in cgs units. The factor fc is the ratio of actual thermal conductivity to 
the classical Spitzer conductivity. The cooling function is a function of temperature T and 
metal abundance Z, and is given by 



A(T, Z) = 2.1 X 10' 
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(6) 
(7) 



in units of ergscm^ s~^. This is an empirical formula derived by fitting to the cooling curves 
calculated by Bohringer & Hensler (1989). We assume that wave injection takes place at radii 
far distant from the cluster center, and thus there is no source term of waves in equation (4). 



The equation for the evolution of shock wave amplitude is given by 



dr 



1 dp 2(7 + l)a^ 2 1 dcs 
p dr CgT r Cs dr 



(8) 



where r is the period of waves, which we assume to be constant (Suzuki 2002). We give the 
period by r = Aq/cso, where Csq is the sound velocity at the average temperature of a cluster 
(Tav), and Aq is the wave length given as a parameter. The second term of the right side 
of equation (8) denotes dissipation at each shock front of the N-waves. We note that the 
sign of the term is the opposite of equation (6) of Suzuki (2002), because we consider waves 
propagating inwards contrary to those in Suzuki (2002). 

For the mass distribution of a cluster, we adopt the NFW profile (Navarro, Prenk, & 
White 1997). The mass profile is written as 



M(r) oc 



In 1 + 



rs(l + r/rs) 



(9) 
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where is the characteristic radius of the cluster. The normahzation can be given by 
-^('^vir) = -^vir, whcre r\i^ and M^i^. are the virial radius and mass of a cluster, respectively. 
We ignore the self-gravity of the ICM. 

2.2. Results 

We show that our model can reproduce observed ICM density and temperature profiles 
of clusters. We choose A1795 and Ser 159-03 clusters to be compared with our model 

predictions. Zakamska & Narayan (2003) showed that thermal conduction alone can explain 
the density and temperature profiles for A1795; /c = 0.2 is enough and other heat sources 
arc not required. On the other hand, the profiles for Scr 159 — 03 cannot be reproduced by 
thermal conduction alone (Zakamska & Narayan 2003). The parameters of the mass profiles 
for the clusters are the same as those adopted by Zakamska & Narayan (2003) and are shown 
in Table 1. The concentration parameter of a cluster, c = ryir/r^ is given by 

1 

c— — 

where Pcrit is the critical density of the universe. We fix the metal abundance profiles. For 
A1795, we assume Z{r) = 0.8exp(-r/170 kpc) Zq (Ettori et al. 2002), and for Ser 159-03, 
we assume Z{r) — 0.51 exp(— r/171 kpc) Zq (Kaastra et al. 2001). 

We carry out the modeling of ICM heating as follows. First, we select values of /c, 
M, and Aq. Then, we set the boundary conditions of the equations (1), (2), (4), and (8) at 
Ti — l kpc, that is, well inside the central cD galaxy. From r — ri,we integrate the equations 

outward and compare the model profiles of ne{ri) and T{ri) with the data. While we fix the 
value of a^(rj), we adjust n^ir,;) and T(rj) to be consistent with the observed profiles. We 
restrict ourselves to a comparison by eye, since neither the data nor the models are reliable 
enough for a detailed fit. If we do not have satisfactory fits, we change the values of /c, 
M, and Aq and repeat the process. We show the values of /c, M, and Aq that we finally 
adopted in Table 2, and briefly summarize how the results depend on the choice of them as 
follows. 

For A1795, we choose /c = 2 x 10^^, because Zakamska & Narayan (2003) have already 
shown that ICM heating only by thermal conduction with /c ~ 0.2 is consistent with the 
observations. In this study, we will show that even when /c is much smaller than 0.2, the 
observed proflles can be reproduced if wave heating is included. However, we found that 
if fc is too small, the obtained temperature profile is too steep to be consistent with the 
observation. For Ser 159-03, we adopt fc — 0.2, which is suggested by Narayan & Medvedev 
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(2001) in a turbulent MHD medium. If we take fc much smaller than this, the model cannot 
reproduce the relatively flat temperature distribution observed in this cluster. 

For mass accretion rates M, we take about 1/10 times the value claimed before the 
Chandra and XMM-Newton era. For A1795, M ~ 500 yr~^ {h = 0.5) was reported 
(Edge, Stewart, & Fabian 1992; Peres et al. 1998). Thus, we adopt M = 5OM0yr-\ which 
is consistent with a recent XMM-Newton observation (< 150 Mq yr~^; Tamura et al. 2001). 
For Scr 159-03, M ~ 300 Mq yr^^ {h = 0.5) was reported (White, Jones, & Forman 1997; 
AUen & Fabian 1997). Thus, we adopt M = 30 M© yr'^ We note that if we assume that 
wave heating is effective and that M is much smaller than the above values, we cannot 
reproduce both density and temperature profiles obtained by X-ray observations; we get too 
high temperature and too low density. 

Typical wave length, Aq, should be comparable to the typical eddy size of turbulence in 
ICM. From numerical simulations, Roettiger, Stone, & Burns (1999) showed that the typical 
eddy size is the core scale of a cluster. Thus, we take Ao = 100 kpc for A1795. For Ser 159- 
03, we use a smaller value of Aq = 70 kpc because of its small mass (Table 1). Smaller Aq 
means a smaller distance that waves propagate before dissipation. 

Among three of the parameters for the boundary conditions at r = {a^, n^, and 
T), we &x aw — 3 to reduce the number of fitting parameters. If we assume much smaller 
a^,, wave heating becomes negligible. On the other hand, if we assume much larger au,, the 
region where the weak shock approximation is invalid {a^ > 1) extends. 

Figure 1 shows the model fits for the two clusters. The boundary conditions are pre- 
sented in Table 2. The temperature T{ri) is especially required to be fine-tuned for the fit. 
We use the Chandra data of A1795 obtained by Ettori et al. (2002) and the XMM-Newton 
data of Ser 159-03 obtained by Kaastra et al. (2001). The XMM-Newton data of A1795 are 
also obtained by Tamura et al. (2001) and they are similar to those obtained by Ettori et al. 

(2002) , although the former is not deprojected contrary to the latter. The good agreement 
between the model and the data suggests that wave heating is a promising candidate of the 
mechanism that solves the cooling flow problem. In Figure 1, densities go to inflnity and 
temperatures go down to zero at r ~ 200 kpc. This suggests that waves injected outside of 
this radius cannot reach the cluster center. 

In Figure 2a, we present the wave velocity amplitude normalized by the sound velocity 
[ayj). As the N-waves propagate into the central regions of the cluster, increases rapidly. 
This is mainly because of the geometrical convergence to the cluster center, whereas the total 
wave luminosity (= energy flux times r^) mostly dissipates through the inward propagation 
in itself. Since > 1 at r < 4 kpc for both A1795 and Ser 159-03, the results may not be 
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quantitatively correct there. We note that the gas velocity for the region of r > 4 kpc is very 
small compared with the sound velocity and thus our ignorance of vdv/dr in equation (2) is 
justified. 

In Figure 2b, the ratio F^fc/Fc is presented. The gaps at ~ 10 kpc reflect Fc < 0. 
Assume that an observer made observations of the model clusters and the temperature 
distributions were exactly measured. If the observer assumed the classical conductivity, the 
heat fiux measured by the observer should be Fd fc {fc < 1) because of the definition of Fc 
(equation [5]). Figure 2b shows that the observer would measure an X-ray emission much 
larger than that predicted by the classical thermal conduction {F^jfc/Fc > 1) if the energy 
swallowed by the black hole at the cluster center is small. Such large X-ray emissions have 
actually been estimated in some clusters (Voigt et al. 2002). Wave heating model can account 
for the observations without the help of heating by AGNs. 



3. One-Dimensional Numerical Simulations 
3.1. Models 

In order to be compared with the results in the previous section, we performed one- 
dimensional numerical simulations. We solve the following equations: 

| + ~(^»=0. (11) 



_„JA(T)-^^1^. (13) 



de 1 r 2 / M ^ d 
_ + _[, vip + e)]^-- 



2 (rj..dT 



where the total energy is defined as e = p/(7 - 1) + and k{T) = KqT^/'^. We Ignore 

the self-gravity of ICM. For numerical simulations, we adopt the cooling function based on 
the detailed calculations by Sutherland & Dopita (1993), 

nlK = [C^{kBT + C^ikBTf + C3]n,ne , (14) 

where rii is the ion number density and the units for ksT are keV. For an average metallicity 
Z = 0.3 Zq the constants in equation (14) are a = -1.7, (3 = 0.5, Ci = 8.6 x 10~^, 
C2 = 5.8 X 10~^, and C3 = 6.4 x 10~^, and we can approximate nine = 0.704(p/mi/)^. 
The units of A are 10~^^ergscm^ (Ruszkowski & Begelman 2002). Note that equation (13) 
does not include the energy dissipation term contrary to equation (4). This is because the 
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energy dissipation at shocks is included automatically in numerical simulations if shocks are 
resolved. 

The hydrodynamic part of the equations is solved by a second-order advection upstream 
splitting method (AUSM) based on Liou & Steffen (1993, see also Wada & Norman 2001). 
We use 500 unequally spaced meshes in the radial coordinate to cover a region with a radius 
of 300 kpc. The inner boundary is set at r = 1 kpc. The innermost mesh has a width of 
~ 15 pc, and the width of the outermost mesh is ~ 3 kpc. The following boundary conditions 
are adopted: 

1. Variables except velocity have zero gradients at the center. 

2. The inner edge is assumed to be a perfectly reflecting point. 

3. The density and pressure at the outermost mesh are equal to specified values. 

Waves are injected at the outermost mesh as 

( 27rCi,oA 

v{t) = a^oCsosm I — — I , (15) 

where a^o is the parameter. 



3.2. Numerical Results 

The mass distribution we assumed is the same as that of A 1795 in Table 1 except for 
Tay. We assume that the ICM is isothermal and in pressure equilibrium at i = 0. For the 
NFW profile (equation [9]), the gas initial density profile is written as 

p(r) = poexp[-5/(r/rJ] , (16) 

where 

m(l) r^kBT.^^ 



m{x) = ln(l + x) — , (18) 

1 + X 

f(x) = l--ln(l + x) (19) 

X 

(Suto, Sasaki, & Makino 1998). The initial ICM temperature is Tav The density and 
pressure of the outer boundary are fixed at the initial values. We finish the calculations 
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when the temperature of some of the meshes becomes zero because we do not treat mass 
dropout from the hot ICM. We define the time as tcooi- If the temperature does not go to 
zero, we finish the calculations at t = 7 Gyr. In this section, we set ne(0) = 0.017 cm~^, 
Tav — 7 keV, and Aq — 100 kpc. Other model parameters are shown in Table 3. 

Figures 3 and 4 show the temperature and density profiles for Models CF and CO, 
respectively. For reference, the data of A1795 are shown (Ettori et al. 2002), although the 
detailed comparison is premature. Contrary to Figure 1, we use a linear scale for the distance 
from the cluster center for both temperature and density profiles to see shock structures. 
Model CF is a genuine cooling flow model; the temperature goes to zero at the cluster center 
at tcooi = 2.7 Gyr. In Model CO, conduction dominates cooling and the solution is stable for 
a long time. The suppression factor of the conductivity, /c = 0.2, in Model CO is the same 
as that of Zakamska & Narayan (2003). The temperature gradient at the cluster center is 
less steep than that obtained by Zakamska & Narayan (2003). This may be because of the 
differences of the boundary conditions or the initial conditions. 

Figures 5 and 6 show the temperature and density profiles for Models Wl and W2, 
respectively. As we predicted in Figure 2, waves are amplified in the central region. A shock 
is seen at r = 13 kpc in Model Wl. Both temperature and density diverge at the cluster 

center because these are one-dimensional simulations and waves focus on the cluster center. 
In the outer region, the shape of waves is complicated. This is because we adopted the 
reflection condition at the inner boundary and the inward and outward waves interact. For 
Model Wl, tcooi = 4.7 Gyr, and for Model W2, tcooi > 7 Gyr. They are much larger than 
^cooi for Model CF. This means the wave heating is effective. 

In Models Ml and M2, we include both thermal conduction and waves. The results of 
Model Ml and Model M2 are almost the same as those of Model Wl and CO, respectively. 
The cooling time-scale tcooi of Model Ml is shorter than that of Model Wl, because shocks 
are weakened by the thermal conduction. In model M2, the conduction dominates the wave 
heating. 

4. Discussion 

Our analytical model and numerical simulations show that acoustic waves are amplified 
at the cluster center and can heat the cluster core. One should note that our assumption, 
that is, the spherical symmetry, could affect the amplification quantitatively. For more 
realistic modeling, we should consider that real clusters are not exactly spherically symmetric. 
However, Pringle (1989) indicated that even if a cluster is not spherically symmetric, the 
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lower temperature and smaller sound velocity at the cluster center should have waves focus on 
the center. In order to study this effect, we need to perform high-resolution multi-dimensional 
numerical simulations (Wada, Fujita, &, Suzuki 2003, in preparation). Since the focusing 
effect depends on the temperature gradient, it may solve the fine-tuning problem of the 
heating of cluster cores. If the cooling dominates heating, the temperature at the cluster 
center decreases and the temperature gradient in the core increases. This strengthens the 
focusing effect and wave heating becomes more efficient. Multi-dimensional simulations 
have another benefit; it is free from the inner boundary conditions that we set in our one- 
dimensional simulations. 

If waves are actually responsible for the heating in cluster cores, weak shocks should be 
observed there in some clusters. In our simulations, waves are amplified at r < 10 — 50 kpc 
(Figures 5 and 6). Thus, if the wave length is large (> 100 kpc), the shocks are not necessarily 
observed in all clusters. However, as pointed out by Churazov et al. (2003), the passing of 
waves may cause gas-sloshing around cluster cores. Observations of fine structures in many 
cluster cores will be useful to understand the heating mechanism there. 

In the central regions of clusters, turbulence generated by interaction of amplified waves 
from various directions could be observed as broadened metal lines by future X-ray obser- 
vatories. Observational studies of optical emission lines revealed that there is warm gas at 
some cluster centers and that the velocity widths range from 100 to 1000 km s^^ (Hu, Cowie, 
& Wang 1985; Johnstone, Fabian, & Nulsen 1987; Heckman ct al. 1989). These warm gas 
may be embodied in and move with the hot turbulent ICM (see Loewenstein & Fabian 1990). 

5. Conclusions 

Through analytical and numerical approaches, we have shown that acoustic waves gener- 
ated by turbulence in the ICM in the outer region of a cluster can effectively heat the central 
part of the cluster. The heat flux by the waves may exceed that by thermal conduction. The 
process presented here is phenomenologically analogous to the collapse of a "tsunami" (a 
seismic sea wave) at seashore owing to the change of depth of the sea. As in tsunamis, even 
if the waves have small amplitude at their origin, they could bring huge damage at a distant 
point, namely the cluster center. Of course, one should note that tsunamis are gravity-driven 
waves and not acousitc waves. In the analj^ical studies, we have obtained time-independent 
solutions and compared the predicted density and temperature profiles with the observed 
ones; they are consistent with each other. In general, we have confirmed the results obtained 
by the analytical studies by one-dimensional numerical simulations. Since we assumed that 
a cluster is spherically symmetric and the assumption leads to artificial focusing of waves. 
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one should take the quantitative results with care. However, it has been indicated that even 
if a cluster is not spherically symmetric, waves are focused by the temperature gradient at 
the cluster center. Thus, it is worthwhile to study the wave heating by multi-dimensional 
analyses. 
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tion, Culture, Sports, Science, and Technology of Japan (Y. F.: 14740 175; K. W.: 15684003). 
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Table 1. Cluster Parameters 



Cluster 




T 

av 




c 




(10^^ Mo) 


(keV) 


(Mpc) 




A1795 


12 


7.5 


0.46 


4.2 


Ser 159-03 


2.6 


2.7 


0.31 


4.7 



Table 2. Model Parameters 



Cluster 


fc 


M 


Ao aw{ri) 


ne{ri) 








(Meyr-i) 


(kpc) 


(cm~^) 


(keV) 


A1795 


2 X 10"^ 


50 


100 3 


0.5 


0.6213 


Ser 159-03 


0.2 


30 


70 3 


0.14 


0.780 



Table 3. Model Parameters 



Model 




fc 


^cooi (Gyr) 


CF 








2.7 


CO 





0.2 




Wl 


0.1 





4.7 


W2 


0.2 







Ml 


0.1 


0.002 


4.3 


M2 


0.1 


0.2 





Note. — No data for tcooi mean 
^cooi > 7 Gyr. 
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/-(kpc) 

Fig. 1. — (a) Modeled temperature and (b) density profiles for A1795 (solid lines) and 
Ser 159-03 (dotted lines). Filled dots and empty triangles are the Chandra data for A1795 
and Ser 159-03, respectively. 




Fig. 2. — (a) Wave amplitudes and (b) the ratio of heat flux by waves to that by thermal 
conduction for A1795 (solid lines) and Ser 159-03 (dotted lines). 
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Fig. 3. — (a) Temperature and (b) density profiles for Model CF. Open circles are the 
Chandra data for A1795 (Ettori et al. 2002). 




Fig. 4. — The same as Figure 3 but for Model CO 
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Fig. 5. — The same as Figure 3 but for Model Wl 




Fig. 6. — The same as Figure 3 but for Model W2 



